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Artifacts arise in the simulations of electrolytes using periodic boundary conditions (PBC). We 
show the origin of these artifacts are the periodic image charges and the constraint of charge neutral¬ 
ity inside the simulation box, both of which are unphysical from the view point of real systems. To 
cure these problems, we introduce a multi-scale Monte Carlo method, where ions inside a spherical 
cavity are simulated explicitly, whilst ions outside are treated implicitly using continuum theory. 
Using the method of Debye charging, we explicitly derive the effective interactions between ions 
inside the cavity, arising due to the fluctuations of ions outside. We find that these effective inter¬ 
actions consist of two types: 1) a constant cavity potential due to the asymmetry of the electrolyte, 
and 2) a reaction potential that depends on the positions of all ions inside. Combining the Grand 
Canonical Monte Carlo (GCMC) with a recently developed fast algorithm based of image charge 
method, we perform a multi-scale Monte Carlo simulation of symmetric electrolytes, and compare it 
with other simulation methods, including PBC-I-GCMC method, as well as large scale Monte Carlo 
simulation. We demonstrate that our multi-scale MC method is capable of capturing the correct 
physics of a large system using a small scale simulation. 

PACS numbers: 82.70.Dd, 82.45.Gj, 05.10.Ln, 87.10.Rt 


I. INTRODUCTION 

Monte Carlo (MC) simulation of Coulomb many body 
systems is a difficult problem mm, mainly due to the 
long range nature of Coulomb interaction. In the most 
naive simulation strategy, one would confine a collection 
of ions inside a box with hard walls, and compute the 
total energy by adding up all pairwise interactions. This 
leads to two obvious difficulties: 1) The artifacts of hard 
wall propagate into the bulk of the system, with the 
characteristic scale set by the Debye length. This ren¬ 
ders large amount of simulation data useless. The situa¬ 
tion becomes worse in the dilute limit. 2) The computa¬ 
tional complexity of each Monte Carlo cycle is of order 
of N'^, with N the system size. This makes simulation 
of large size systems practically impossible. Obviously, 
both difficulties have their roots in the long range nature 
of Coulomb interaction. 

Adoption of periodic boundary conditions (PBC) re¬ 
stores the translational symmetry and hence eliminates 
all the boundary effects. Furthermore, one can take ad¬ 
vantage of periodicity and compute the long range part 
of Coulomb energy in Fourier space, using the method 
of Ewald summation [3]. (For a pedagogical introduc¬ 
tion, see [ 2 ].) This reduces the computational complexity 
from to which is still not fast enough for large 

scale simulations. For molecular dynamics (MD) simula¬ 
tions, one can use particle-mesh method and fast Fourier 
transform n m to reduce to complexity further down 
to TV log A^. Unfortunately, these techniques are not ap¬ 
plicable for Monte Carlo simulations. This is probably 
the main reason why MD has become the more popu¬ 
lar method for charged systems, even when one is only 


insterested in static properties. 

By periodically replicating the system, one introduces 
an infinite array of image charges for each charge be¬ 
longing to the system under study. These images are 
unphysical and break the rotational symmetry. Further¬ 
more, the total charge inside the simulation box must 
vanish, for otherwise the summation over images would 
diverge. This leads to an unphysical constraint of charge 
neutrality inside a finite volume. In Fig. 1(a) we show 
the average charge density around a fixed ion, obtained 
using different simulation methods. Whilst large scale 
MC simulation (STD) gives a clean form of screened 
Coulomb, small scale simulation (with system size three 
Debye lengths) using PBC (PBC-I-GCMC) gives substan¬ 
tially different results. The deviation, which grows with 
the distance to the test ion, is caused by the unphysi¬ 
cal image charges. Similarly, internal energies calculated 
using PBC also deviate substantially from the standard 
results, shown in Fig. |l(b)| These artifacts are quan¬ 
titatively less important for larger systems, and become 
invisible for infinite systems. In principle, one can always 
simulate sufficiently large systems so that these artifacts 
become insignificant, or correct these artifacts for every 
microscopic configurations during the simulation. This 
however generically leads to waste of computational re¬ 
sources, or slow-down of simulation processes. 

As an alternative to PBC, one may use a linear con¬ 
tinuum theory to model the influence of the subsystem 
outside the simulation domain, whilst particles inside the 
domain are simulated explicitly. For a dipolar system, 
the subsystem outside responds to the dipoles inside the 
cavity and exerts a field on the latter, which are con¬ 
ventionally called the reaction field. The term reaetion 
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FIG. 1: (a) Screening charge density around a fixed ion, as a function of distance to the source. Our multi-scale MC method 
(dark triangle, RP-I-GCMC) agrees the standard results (red dot, STD) to a high precision, indicating that it correctly captures 
the physics of an infinite size system. By contrast, simulation using PBG (pink diamond, PBG-I-GGMG) gives substantially 
different results. The straight line is prediction by LPB in open space, Eq. (551, whereas the dashed curve is the prediction by 
LPB with PBG, Eq. System parameters in these simulations are listed in Sec. IV B (b) The average energy per ion, in the 
unit of feeT, as a function of ion density. Again our multi-scale MG method agrees with the standard results to a high precision, 
whilst simulation using PBG yields largest errors among all methods. The solid line (DH) is the prediction of Debye-Hiickel 
theory in open space (the first term in RHS of Eq. and the dashed line, Debye-Hiickel theory with PBG (Eq. (|10[)). 



FIG. 2: Schematic illustration of the multi-scale reaction 
potential grand canonical Monte Garlo. If this simulation 
method work well, the boundary of simulation cavity must 
behaves as a virtual one, i.e., the simulated system must be¬ 
have as ions inside a virtual spherical domain in an infinitely 
large system. 

field boundary condition is named after Onsager , who 
first used this method to calculate the dielectric constant 
of dipolar fluids. Similar multi-scale ideas, i.e. treating 
dipolar molecules in near and far fields using different 
methods, however can be dated back to Clausius [7] and 
Mossotti [8] . For a discussion, see the monographs on di¬ 
electrics by Frdhlich [3], and by Battcher m- Born m 
also used a similar idea to calculate the solvation (free) 


energy of ions in solvent. We note that in the study of 
ionic systems, the term reaction potential would be more 
appropriate than reaction field, since it is the electrostatic 
potential, rather than field, that couples to the charges. 
The same method can be generalized to molecules with 
arbitrary charge distributions. 

Simulation methods using multi-scale strategy have 
been studied by many authors [T^2()j . These methods 
were designed for study of either dipolar systems or ionic 
systems. The simulation cavity can be either fixed and 
common for all simulated particles, or can be moving 
and individual for each particle. The subsystem outside 
the cavity is always assumed to obey Poisson equation 
(for dipolar systems) or the linearized Poisson-Boltzmann 
equation (for ionic systems). Finally, it appears that in 
most of the previous works, the reaction-field modeling 
is used in molecular dynamics simulations (MD), even 
though in principle, it is applicable for Monte Carlo sim¬ 
ulations as well. 

To implement the multi-scale simulation, one must 
solve the continuum theory for each simulation step. The 
classical Kirkwood expansion [5T], though straightfor¬ 
ward, is not efficient because of its slow convergence. The 
image charge method, which was first discovered by Neu¬ 
mann |22j for polar systems, and have been generalized 
to ionic systems , transforms the Kirkwood expan¬ 

sion into a line integral, and therefore can be efficiently 
computed. Recently, we improved the image method [26] 
(using Mellin transform technique) so that it is applicable 
for Poisson-Boltzmann theory. 

Regardless of the rich history, there are still a few con- 
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tributions that we can make to improve this multi-scale 
simulation strategy. Firstly, in order to mimic an infinite 
system using a finite and simulation domain (which is the 
essential task of almost all computer simulations), the to¬ 
tal charge inside the simulation cavity must be allowed 
to fluctuate. This means that the numbers of positive 
and negative ions must be allowed to fluctuate indepen¬ 
dently. A “grand canonical Ewald” simulation where the 
ions are inserted and deleted in neutral pairs can not 
capture the physics of charge fluctuations. In this work, 
we shall discuss in detail the artifacts due to suppression 
of charge fluctuations. Secondly, all previous works as¬ 
sume without proof that the subsystem outside the cavity 
is described by the linearized Poisson-Boltzmann theory 
(LPB) for ionic systems. We shall derive the correct lin¬ 
earized continuum theory and show that it reduces to the 
naive LPB only in the dilute limit. For non-dilute elec¬ 
trolytes, the continuum theory may be considerably com¬ 
plicated. Furthermore, for asymmetric electrolytes, there 
is also a constant cavity potential that acts on all ions in¬ 
side the cavity. Thirdly, we shall combine the multi-scale 
GCMC with a recently developed image charge method 
[25] to simulate symmetric electrolytes. We find that 
multi-scale simulation of a small system (with linear size 
of only three Debye length) is capable of capturing the 
physics of an infinite system. For example, as shown in 
Fig. |l(a)| and Fig. |l(b)[ both the screening charge density 
around a fixed ion and the average internal energy per ion 
computed using our multi-scale GCMC (RP-b GCMC) 
agrees remarkably well with the standard results (STD). 
By contrast, to achieve similar precision in simulations 
using PBC, the system would have to be at least as large 
as ten times Debye length. 

The remaining of this paper is organized as follows. 
In Sec. jn] we analytically discuss LPB with PBC, and 
demonstrate two sources of artifacts of PBC, i.e., periodic 
image charges, and constraint of charge neutrality inside 
each unit cell. In Sec. |III[ we start from the basic prin¬ 
ciple of statistical mechanics, and systematically derive 
the multi-scale MC simulation strategy. This derivation 
made it clear that one should always use grand canonical 
ensemble for the subsystem inside the cavity, and that, 
for asymmetric electrolytes, there is a constant cavity po¬ 
tential acting on every ions inside. In Sec. |IV| we outline 
the numerical implementation of the multi-scale strategy 
using grand canonical MC simulation, as well as that of 
other simulation methods that are used for comparative 
studies. In Sec. jVj we simulate the symmetric primitive 
model, and compare the simulation results using different 
methods. We demonstrate how our multi-scale method 
captures correctly the correlation effects between ions, 
whilst the Ewald summation methods fail to do so. Fi¬ 
nally in Sec. |VI| we draw the conclusion remarks. 

To simulate large systems, one need to use fast algo¬ 
rithms based on multipole expansion, such as the oct- 
tree algorithm [27], to expedite the computation of elec¬ 
trostatic energy. A Monte Carlo simulation strategy 
that combine multi-scale modeling, grand canonical en¬ 


semble, image charge techniques, and oct-tree algorithm 
speeding-up can achieve a computational complexity of 
order of A^logfV for one Monte Carlo cycle, and is free of 
artifacts due to PBC. It therefore constitutes a competi¬ 
tive alternative to molecular dynamics simulation meth¬ 
ods for electrolyte systems. This will be discussed in a 
future work. 


II. ARTIFACTS DUE TO PERIODIC 
BOUNDARY CONDITIONS 

Consider N positive ions q and N negative ions —q 
inside the simulation cube, which is centered at the origin 
and has length L, and volume V = L^. To simplify the 
analysis, we shall only focus on symmetric electrolytes in 
the dilute regime, where linearized Poisson-Boltzmann 
theory is applicable. We shall for the moment assume 
that all ions that point-like, and consider the case of finite 
ion sizes towards the end of this section. We periodically 
replicate the system, so that both the charge distribution 
and the electrostatic potential become periodic functions. 

We fix a positive ion q at the origin, and let all other 
ions fluctuate according to the equilibrium Gibbs mea¬ 
sure. We shall calculate two quantities: 1) the screening 
charge density Pq{r) around the test ion, and 2) the cor¬ 
relation potential ip, i.e., the mean potential acting on 
the test ion due to all other ions. These two quantities 
are calculated in the classical Debye and Hiickel theory 
(DH) [30] for an infinite electrolyte. What we are doing 
here is to generalize DH to a finite system with PBC. 
The total internal energy of the system is related to ip 
via 



Now according to the linearized Poisson-Boltzmann 
theory (LPB), the average number densities of mobile 
positive and negative ions are related to the mean poten¬ 
tial (p{v) via the Boltzmann factor: 

c-e(i') = « c+- ^gc+(/)(r), (2a) 

c_(r) = « c_-b/3gc-0(r), (2b) 

where we have linearized the exponentials in the second 
steps. For a physical system with infinite size, one usu¬ 
ally requires that the mean potential (p{v) vanishes in the 
infinity. Then according to Eqs. Q, c± are the number 
densities of positive/negative ions in the infinity, i.e., the 
bulk ion densities [55] . 

This argument is no longer applicable if periodic 
boundary condition is imposed, since the mean potential 
(piy) does not vanish as infinity. For periodic systems, the 
most convenient choice is that the mean potential (p{v) 
vanishes identically when integrated over the unit cell: 
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This is equivalent to the “conductor boundary condition” 
that are popularly used in Ewald summation method, 
which requires that the Fourier transform of electrostatic 
potential does not have a k = 0 component. Now if we 
integrate Eqs. (|^ over the unit cell with volume V, the 
LHS are the total numbers of mobile positive/negative 
ions in the unit cell, which are JV — 1, and JV, respectively, 
whereas the RHS are c±V respectively. This allows us to 
determine the parameters c± as 


test ion itself from Eq. Q, and take the local limit r —>■ 0: 


V' = 


lim 

r—i-O 



g 

47re|r| 


QK qe '‘I'l 

47re 47re|l| 

MO 


g 

e/i^V' 


(9) 


c+ 


JV-1 _ N 
V ’ ~ V 


( 4 ) 


Let us now invoke the Poisson equation, use Eqs. ([^, 
and expand to the first order in ij). We obtain 


- A(j){r) = -p,(r) =-g Vd(r-1) + g(c+(r) - c_(r)) 
e e “ 

= (5) 

1 


where 1 = {L^rix, LyUy, L^riz) are lattice vectors, n^,, Uj,, 
riz are integers, and k is the inverse Debye length: 


K = A/e-i/3g2(2Ar- 1)/E. (6) 


The first term in RHS of Eq. ([^ is due to the test ion and 
its periodic images, whilst the second term corresponds to 
a uniform and negative charge density, arising due to the 
fact there is one more negative mobile ion than positive 
mobile ions, see Eq. Q. This in turns originates from 
the artificial constraint of charge neutrality inside each 
primitive cell. Integrating both sides over a unit cell, 
one easily see Eq. ^ is indeed satisfied. Note that inverse 
Debye length k in Eq. Q is also slightly different from 
its usual form in free space. 

Eq. (§ can be easily solved subject to PBC: 


1 


gg-ft|r-l| 

47re|r — 1| 


g 


( 7 ) 


which can be easily shown to satisfy the condition Eq. (|^ . 
Taking the Laplacian of this equation, we find the total 
average charge density: 


The first term is just the correlation potential in free 
space as predicted by Debye-Hiickel theory, and the other 
two terms constitute the artifacts due to periodic bound¬ 
ary conditions, both of which vanish as the system size 
L goes to infinity (with the Debye length fixed). 

If the ions are not point-like, but are hard spheres with 
diameter d (the primitive model), the above results need 
to be modified properly. The resulting theory is consid¬ 
erably more complicated. If the ion density is not very 
high, however, we may consider corrections due to hard¬ 
cores as small perturbations, same as the artifacts due to 
PBC. Then all we have to do is to correct the first term 
in Eq. ([^, i.e., replace it by the correlation potential of 
a hard sphere ion, which was worked out by Levin and 
Fisher some time ago [551 US]- This leads to 


V- = - 


q K 


47re(I ■ 


Ka) 


E 

MO 


qe 


-k|1| 


47re|l| 


CK^V 


( 10 ) 


This result is plotted in Fig. |I(b)| (PBC-I-DH), which 
agrees reasonably well with the MC simulations with 
PBC if the ion density is low. The agreement becomes in¬ 
creasingly worse as the density increases, indicating that 
the approximation underlying Eq. (101 becomes increas¬ 
ingly inaccurate. 

In Fig. [2 we have chosen the size of simulation box to 
be approximately three Debye length, in order to high¬ 
light the artifacts due to PBC. In many simulations, 
the size of simulation box is chosen to be much larger 
than the Debye length. In this case, the second term in 
Eq. (10) is exponentially small comparing with the third 


term. The artifacts is therefore completely due to the last 
term. The relative error (in the dilute regime ha —> 0) 
due to PBC is the given by 


“E 


qhf e '^1'’ ’ 
47r|r — 1| 


( 8 ) 


dip 47r 

Ip [RLY 


( 11 ) 


where the first term is due to the test ion and its images, 
whilst the second term is sum of the screening charge 
densities around these ions. Eq. (|^ (ave raged over all 
orientations of r) is plotted in Fig. 1(a) as PBC-I-DH, 


which agrees remarkably well with MC simulations with 
PBC (PBC-bCCMC). 

The correlation potential ip acting on the fixed ion can 
be easily obtained by subtracting off the part due to the 


To make this relative error less than 1%, one would have 
to choose the size of simulation box to be larger than 
ten times Debye length. This may be very difficult to 
achieve in many simulations. By strong contrast, using 
the multi-scale GCMC simulation method, we can faith¬ 
fully simulate an infinite system using a simulation cavity 
with radius only three times Debye length. 
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III. MULTI-SCALE MODELING 
A. Effective Interaction 


where Y) is the potential acting on the charge qi at 
Xi due to all ions in B, in the specific micro-configuration 
Y = {yi,...,yM}: 


We shall study the primitive model of electrolytes, 
where the solvent (water) is modeled implicitly as a di¬ 
electric medium with relative permittivity e/eo = 80, 
whilst ions are modeled as hard spheres with a point 
charge at the center. Furthermore, the permittivity of the 
hard spheres is assumed to be the same as that of the sol¬ 
vent. We introduce a spherical simulation cavity, which 
is already schematically illustrated in Fig.[^ The subsys¬ 
tem A inside the cavity consists of N charges {qi ,..., q^} 
with positions denoted by X = {xi,..., xjv}, whereas 
the subsystem B outside the cavity consists of M charges 
{Ox, ..., Om} with positions denoted Y = {yi,..., yM} 
respectively. 

The total Hamiltonian of the overall system consists of 
three parts: 

H = + + (12) 

where and are the Hamiltonian of isolated sub¬ 
systems A and B: 




E 


i<3 


1 qiOj 
47re |xi - Xjl 


-I- w(Xi - Xj) 


E 


OL<^ 


1 QgOp 
47re |ya - yp\ 


+ w(ya 



(13a) 

(13b) 


The first term in each bracket is the long-range Coulomb 
interaction and the second term w(x — y) is the short 
range interaction, which is assumed to be independent of 
the ion species . For the primitive model, w(x — y) is 
just the pairwise hardcore repulsion: 


M 

<p(x„Y) = 

a—1 


1 0a 

47re |xi - yal’ 


(17) 


The grand canonical partition function of the overall 
system can be written as 

S = Tr^TrBe-^(^^+^"+"^'’), (18) 

where the traces Tr^ and Trg represent integrations of all 
the variables as well as summation over particle numbers: 



Here Ng and Mg are the numbers of mobile particles of 
s-th species in subsystems A an B respectively, whereas 
/is are their chemical potentials, and Kg are some micro¬ 
scopic length scales. 

We integrate out the variables in subsystem B, and 
obtain an effective theory for the subsystem A'. 

S = Tr^e-^^^(TrBe-''K+^^")) 

= Tr^ = Tr^ ^ (21a) 

where 6H^ represents the additional effective interac¬ 
tions of subsystem A as mediated by subsystem B: 


;(x-y) = 


0 , 


|x-y| > d; 


oo, |x-y| < d. 


(14) 


The last part in the RHS of Eq. (12), represents 


the interaction between subsystems A and B: 
N M r . 


H 


AB 


= EE 

i—l oc—1 


47re |xi - yo 


+ w{x, - y„) 


(15) 


We note that the short range part of in Eq. (15) 


becomes nonzero only if the distance between Xi and y^ 
is smaller than d, which can happen only if both particles 
are very close to the cavity boundary, hence it has no 
influence on the ions that are at least one ion-diameter 
away from the boundary. We shall therefore ignore this 
short range interaction between A and B. Accordingly, 
we shall exclude a thin spherical shell near the cavity 
boundary for collection of simulation data, id^® can 
then be written into the following form: 


N 




( 16 ) 


e-^'^"=TrBe 




(21b) 


Therefore we can study the subsystem A only with an 
effective Hamiltonian 

= + SH^ = H^- (22) 

and the influences of the subsystem B on A are auto¬ 
matically taken care of. We will find that can be 
approximately calculated and the result assumes a simple 
and physically transparent form. 

We emphasize that Eqs. (21a I and (18) are grand 
canonical partition functions. Hence the ion numbers 
for each species in A and B are all stochastic variables. 
Now it is well known that for large systems (such as sub¬ 
system B), all ensembles are equivalent to each other. 
But for small systems (such as A), different ensembles 
are inequivalent. It is therefore important to use grand 
canonical ensemble, rather than canonical ensemble for 
A in our multi-scale modeling. We shall demonstrate this 
point by comparing simulations results using GCMC and 
CMC in Sec. El 
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B. Method of Debye Charging 

To calculate the effective interaction 5H^, we shall use 
the Debye charging method [30]. Let us scale charges of 
all ions qi in subsystem ^ by a factor A: 






SO that Eq. (16) becomes 

N 




AB 


i=l 


Eq. (21b) then becomes: 


(23) 


(24) 


(25) 


Taking derivative of Eq. ( |25| ) with respect to A, we find 

(26) 


N 




i=l 


where the A-dependent average {0)x is defined as 
TrBOe-/5(^®+^«^^) 


(0)a = 




(27) 


The physical significance of ((/7(xi,Y))^ is therefore the 
average potential at (inside the cavity, occupied by an 
ion Xqi) due to all ions in B. Note that the average is 
over the statistical fluctuations of all ions outside, with 
the locations of all ions Xqi inside fixed. To calculate 
this average, we treat A as a small parameter and use the 
linear response theory. That is, we expand Eq. (27) to 
the first order of A: 

(0)a « (O)o - A/3 - (O)o ■ 

+ 0(A2) (28) 


Note that the average (O)o in RHS is defined by Eq. (27) 


with A set to zero, i.e., with the interaction between two 
subsystems switched off: 


(O)o = 


Tr« ' 


(29) 


Let us define a kernel x(r, r') in terms of the connected 
correlation function of :^(r, Y) as follows: 

X(r,r')= - /3[(v3(r,Y)(/j(r', Y))o 

- (<7>(i’,Y))o ((/j(r', Y))o]. 


Setting O = ip{-x.i,Y) in Eq. ( |M| ) and using Eqs. (29), 
(16), we obtain: 


{^p{:>^i,Y))x = ((p(x„Y))o + A^gj-x(xi,Xj)- (30) 


The first term in the RHS is the average potential at x^ in 
the absence of any charges in A, generated by an isolated 
subsystem B. We shall call this term the cavity potential 
‘f’cav(r), following the terminology of Onsager |31|. Evi¬ 
dently it satisfies the Laplace equation inside the cavity, 
and is invariant under arbitrary rotation. It then follows 
that $(.av(r) must be a constant in the whole cavity: 


$cav = (¥>(r,Y))^ 


(31) 


Furthermore, for symmetric electrolytes, the cavity po¬ 
tential vanishes identically due to charge inversion sym¬ 
metry. For asymmetric electrolytes, $cav is generally 
non-vanishing. This term has been missed by all previous 
works on reaction-field modeling of charged systems. 

The second term in the RHS of Eq. (30) is linear in A. 


It therefore arises due to the linear reaction of subsys¬ 
tem B to sources charges {A gi,..., AgAr} in Al. We shall 
therefore call it the reaction potential |35j . 

Finally, combining Eqs. (31) and (30), and substituting 


them back into Eq. (261, integrating the latter over the 


charging parameter A from 0 to 1, we obtain the effective 
interaction of A mediated by B: 


N 




(32) 


bi=i 


The correlation function x(r,r') is intimately related 
to the electrostatic Green’s function of the subsystem B. 
To see this, let us insert a test ion with charge g at r' 
inside the otherwise empty cavity, and calculate the av¬ 
erage total potential at r, which may be inside or outside 
the cavity. This potential is the linear superposition of a 
part due to the source charge g, and another part due to 
all ions in subsystem B. The later has the form given by 
Eq. (30), to the first order in g (with A set to unity, of 
coursejT Hence the total mean potential at r is given by 


dlcav T g 


1 




0(g2). 


(33) 


The sum inside the bracket describes the linear response 
of mean potential at r due to a unit point charge at r', 
and therefore is the electrostatic Green’s function: 


G(r,r') = 


47re|r — r' 


Tj +X(i',r')- 


(34) 


Substituting this back into Eqs. (22) and (32), we find 


the effective Hamiltonian 


N 


1 


i<j 


E 


g*$cav + Ag2^(Xi,Xj) 


(35) 


where the first term is the effective interaction between 
different ions, and the second term is the self-energy of 
ions mediated by subsystem B. 
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C. Linearized Poisson-Boltmann Theory 


The Green’s function G(r,r') encodes the linear re¬ 
sponse properties of the subsystem B, that is an elec¬ 
trolyte with a spherical void (the empty cavity). More 
precisely, it is the incremental mean potential at the field 
point r due to the insertion of a unit test charge at the 
source point r'. If r is inside the cavity, G(r, r') obviously 
satisfies the Poisson equation: 

— ArG(r, r') =-(5(r — r'), r & A. (36a) 

If r is outside the cavity, it satisfies a linear integro- 
differential equation: 


— ArG(r, r') -I- J a(r, r")G(r", r')dr" = 0, r £ B. 

(36b) 

The kernel Q!(r,r") is generally unknown. At the level 
of first order perturbation, however, we may approxi¬ 
mate Q:(r,r") by the kernel that corresponds to a uni¬ 
form electrolyte, which only depends on the difference 
r — r", because of the translational symmetry. The 
Fourier transform of this kernel was calculated analyt¬ 
ically recently, both for symmetric and asymmetric elec¬ 
trolytes [32] • The results are still quite complicated. 

In this work, we shall only study symmetric elec¬ 
trolytes in the dilute regime, where LPB is applicable 
in the bulk. This entails two essential simplifications for 


the effective Hamiltonian Eq. (35). Firstly the cavity po¬ 


tential vanishes identically. Secondly, the kernel a(r, r') 
has the simple form 6{r — r'), where k is the inverse 
(bare) Debye length, given by 


= e-S 




(36c) 


with c the average ion density. Consequently, Eq. (36b) 


reduces to the linearized Poisson-Boltzmann equation: 

- ArG(r,r')-f K2G(r,r') =0, r G H. (36d) 


We shall deal with the cases of dense electrolytes and 
asymmetric electrolytes in a future work. 

We also need to determine the boundary conditions 
satisfied by the Green’s function. At infinity, it clearly 
satisfies the free boundary condition: 

lim G(r, r') —)■ 0. (37a) 

r—>-oo 

On the cavity interface, r = R, G(r,r') and its normal 
derivative are continuous: 


lim G(r,r') = lim G(r,r'), (37b) 

r—^R’' 1 — yR+ 


lim 


9G(r, r') 
dr 


lim 

1 — yR+ 


dG{r, r') 
dr ’ 


(37c) 


where r —> mean that the field point r approaches 

the interface from outside/inside, respectively. 


Our method therefore works as follows. Firstly we find 
the Green’s function by solving Eqs. (36a), ( 36d ), subject 
to boundary conditions (37), then use Eq. (34) to find the 
reaction potential x(r, r^Tand finally use Eq. (351 to find 
the effective Hamiltonian. 


IV. NUMERICAL IMPLEMENTATION 

In this section, we discuss the numerical implemen¬ 
tation of our multi-scale MG method, as well as other 
simulation methods that we use for comparison. 


A. Computation of the Green’s Function 


In one of our previous works [23], we discus sed a n ef¬ 
ficient image charge method for solving Eqs. (36a) and 
(36d), subjected to boundary conditions Eqs. (37). Here 


we briefly summarize the results in order to make this 
work self-contained. Let r' and r be, respectively, the 
source point and the field point. We obtain the reaction 
potential x(r, r') in the form of Kirkwood series [2Ij : 


X(r,r') = 


I 


oo / / \ n 

rr ^ 


AttcR 


n—0 


M„(m)P„ (cos 0), (38) 


where 9 is the angle between r and r', Pn{x) the Legendre 
polynomials, R the cavity radius, u = kR, and 


M„(m) = 


(n -I- l)kn{u) + uk'^{u) 


nkn{u) — uk'^{u) ’ 

with kn{u) the modified spherical Hankel functions: 

{n + iy. 1 


kn{u) = 


ire 


2u f^^l\{n-l)\(2uy 


E 


(39) 


(40) 


Eq. (38) has a useful integral representation: 

X{t) 


, „ I Xit) , 

= pry'*. 


(41) 


where the vector t is parallel to r' and can be written as 


t = ti - ] = tr\ 


Its magnitude t runs from vk = R? /r' to inhnity. Eq. (41) 


has an intuitive physical interpretation: it is the poten¬ 
tial at r due to a fictitious straight-line of image charges 
which starts at the Kelvin point tk = R?r' jr'"^^ and ex¬ 
tends all the way to infinity. As shown in reference [26], 
the line charge density A(t) [36] is related to the Mellin 
transform of M„(rt) (as a function of n), and is an oscil¬ 
latory function of t. 



































Line image 
charge:A(x) 


FIG. 3: Schematics. A source charge inside the cavity (red 
point), and the reaction potential x)!", r') generated by a line 
of image charges (blue straight line) as given by Eq. ( |41| |. The 
blue color inside the cavity is the density plot of screening 
charge density, obtained directly from simulation data. The 
dashed circles are the contour lines as predicted by LPB. The 
fact that these two plots agree with each other demonstrates 
the physics of an infinite system is faithfully reproduced in 
our multi-scale simulation using a small cavity. 


size is not much larger than the Debye length, and it is 
a nontrivial matter to cancel the influences of boundary 
and restore the physics of an infinite system. All ions 
have hardcore diameter d = 7.5A. 

We simulate the subsystem A with an effective Hamil¬ 
tonian Eq. (35), with the reaction potential computed us¬ 
ing Eq. (42). As emphasized in the preceding section, we 
must use grand canonical Monte Carlo to simulate this 
system. Use of canonical Monte Carlo would suppress 
fluctuations of total charge inside the cavity and there¬ 
fore leads to substantial errors. We will demonstrate this 
point in great detail below. 

The probability density function of a microstate with 
N ions {qi ,..., q^} at X = {xi,...,xv} is given by the 
grand canonical Gibbs distribution: 


N 


Pn 


(X)[]d3x, = ^ 




n 


iV'A: 


3Ns 


e 


-0H 


A 

eff 


N 


i 


^,(44) 

Note that the grand canonical partition function S is di¬ 
mensionless, whilst the dimension of p(X) is . To 
perform numerical simulation, however, it is mandatory 
to deal with discrete probabilities. Therefore we dis¬ 
cretize the simulation domain, with elementary length 
unit |dx|. The probability of a discrete microstate is then 


To further reduce the computational cost, we truncate 
the line integral and discretize the line image using Gaus¬ 
sian quadrature. can then be expressed by M 

point-like image charges and a few correction terms: 


M 


m—0 ' 


^XiPiicosO) (42) 


where M image point charges are located at r^, with 
magnitude qm. The correction terms are due to the trun¬ 
cation of the line integral. Numerical tests showed that 
4 images and L = 1 corrections provide result with error 
less than 1% in computing the self-energy of a charge. 
All relevant details of qm, I'm and xi can be found in ref¬ 
erence |26] . For not too small r, Eq. (421 is much more 


efficient than the Kirkwood expansion, Eq. (38). 


If the source point approach the center of the sphere, 
r R, the Kirkwood expansion converges rapidly. We 
can use leading term of the Kirkwood expansion to ap¬ 
proximate the Green’s function: 


X(r,r' = 0) = - 


47re(l -I- kR) 


(43) 


B. Grand canonical Monte Carlo 

The system parameters used in our multi-scale simu¬ 
lation are as follows. The radius of simulation cavity is 
fixed to be lOOA. The total ion number (including both 
species ) varies from 15 to 40. The corresponding range 
of Debye length is between 24 — 50A. Hence the system 


^iv(X) 


1 


e 


-0H 


A 

eff ^ 


(45) 


The Grand Canonical Monte Carlo simulation consists 
of three steps: displacement, insertion and removal 

(1) Displacement. We select an ion randomly, and let 
its charge and position be qt and x^. The part of total 
energy of the system that depends on is given by 




|^g*gjG(xi, Xj) -I- w(xi - x^) 

3¥=i 


fA(x„ 


(46) 


Then we choose a destination x- for qi according to a flat 
probability density function defined in a cubic cell with 
a given size and centered at x^, and let the new energy 
be E[. The change of the total energy of the system is 

AE = E[-Ei. (47) 

The displacement is accepted with probability 

acc(X->X') =min{l,e"'^^-®} (48) 

Evidently, if any hardcore constraint is violated in the 
destination state, AE = oo, and the displacement is re¬ 
jected with probability one. If the center of an ion across 
the boundary of simulation domain in an attempted 
move, it is also rejected. In another word, the bound¬ 
ary behaves as a hard wall to the center-of-mass of all 
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ions. It is easy to see that this transfer probability satis¬ 
fies the detailed balance with respect to the equilibrium 
Gibbs distribution Eq. (45): 

7 r(X) acc (X ^ X') = 7 r(X') acc (X' ^ X). (49) 


(2) Insertion. To mimic the fluctuations of particle 
numbers in an open system, we need to insert and remove 
particles. In a micro state with Ng particles of species 
s, we first randomly select a position x in the cavity ac¬ 
cording to a uniform probability distribution. The latter 
can be generated by introducing three independent ran¬ 
dom variables ryi, 772 , 773 that are uniformly distributed in 
the interval [ 0 , 1 ], and express the spherical coordinates 
of the insertion point x as 

r = rj\^^R, 6* = cos“^(l - 2772 ), (j) = 2TTr]3. (50) 

It is straightforward to verify that the Jacobian 
9 (a 7 ,7/, z)/ 9 (? 7 i, 772,773) is a constant independent of 771, 
V 2 , where x,y,z are the Cartesian coordinates. We 
now choose a species s randomly and insert an ion at 
X with probability 


acc {Ng Ng + 1) = min < 1 , 


Mnl-Ei) 


Ng 


1 


(51) 


where Ei is the energy change due to the insertion op¬ 
eration, given by Eq. (46), whilst y* are the effective 
chemical potential of the sth species , given by 

V 


fj,* = yg+ feeTlog 


A?- 


(52) 


C. Comparison with Other Simuation Methods 


Three other simulation methods are used to compare 
with our multi-scale reaction potential GCMC simula¬ 
tion method (RP-I-GCMC). To demonstrate the artifacts 
due to periodic boundary conditions, we use the popu¬ 
lar Ewald summation GCMC method (PBC-I-GCMC) to 
conduct a small scale simulation. To demonstrate the ar¬ 
tifacts due to suppression of charge fluctuations, we apply 
a small scale simulation on reaction potential model using 
canonical Monte Carlo (RP-I-CMC). Finally, we conduct 
a large scale simulations with system sizes at least ten 
Debye lengths to obtain standard results (STD), with 
which all other simulations are compared. To make the 
comparison meaningful, we adjust parameters such that 
the ion densities in the bulk are equal in all simulations. 
Furthermore, all small scale simulations have equal size 
of simulation domains. 


(1) Grand canonical Ewald summation MC with Peri¬ 
odic Boundary Conditions (PBC-hGCMC). We choose 
cubic simulation domain with volume equal to 47 ri?^/ 3 . 
To avoid divergence when summing over periodic image 
charges, the total charge in the simulation box must van¬ 
ish in every micro state. Furthermore, one still need to 
specify some “boundary conditions” at infinity. We shall 
use the popular “conductor boundary condition”, where 
the average potential inside the simulation box vanishes, 
and there is no dipolar term in the total free energy. For 
details, see the textbook by Frenkel and Smit [2]. This 
choice is consistent with the condition Eq. that we 
use in the analysis of LPB with PBC. 


(3) Removal. In a micro state with Ng particles of 
species s, we randomly choose a particle with (a ran¬ 
domly chosen) species s, and remove it with probability 

acc {Ng Ng — 1) = min |l, Ng (53) 


where —Ei is the energy change of the removal operation, 
with Ei again given by Eq. (46). 


It is straightforward to verify that the probabilities 
of insertion and removal satisfy detailed balance with 
respect to the equilibrium grand canonical distribution 
Gibbs distribution Eq. (45): 


TTjv I ) acc {Ng —>■ Ng 1 ) 


— ’’’Af-i-i 1 ^ ^ ) 3'CC {Ng 1 ^ Ng). (54) 


The factor jJxp/F in the LHS is the probability of choos¬ 
ing one particular point inside the (already discretized) 
simulation domain, whereas the factor l/{Ng -j- 1) in the 
RHS is the probability of choosing a particular particle 
with species s. 


The total ion number however can fluctuate accord¬ 
ing to the standard grand canonical ensemble theory. In 
accordance with these constraints, ions are inserted and 
deleted at random as pairs [T]. 

(2) Multi-scale Reaction Potential Canonical MC sim¬ 
ulation (RP -h CMC). This simulation is similar to 
RP-I-GCMC, with the only difference that the ion num¬ 
ber of each species is kept constant. 


(3) Large Scale Canonical MC simulation with Hard Wall 
Boundary Conditions (STD). We conduct a large scale 
canonical simulation using a spherical simulation domain 
with radius at least ten times of Debye length, and with 
hard wall boundary condition. To void the artifacts due 
to boundary effects, only ions that are more than five De¬ 
bye lengths away from the boundary are used for data col¬ 
lection. The simulated system typically contains about 
2500 to 5000 ions. To speed up the computation of total 
electrostatic energy, we use the recently developed oct- 
tree algorithm [27]. We calculate various physical quan¬ 
tities via this method and use them as standard results 
for all later comparisons. 
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FIG. 4: Ion density profile inside the simulation cavity, as a 
function of radius. The density profile of RP+GCMC is al¬ 
most flat inside the cavity, except in a boundary layer of ap¬ 
proximately one ion diameter thick. By contrast. The profile 
of RP-I-GMG shows a much larger variation inside the cavity. 
This is caused by the artificial constraint of charge neutrality 
inside the cavity. The standard results are obtained using a 
large scale simulation, and therefore are completely flat in the 
plotted region. 


V. RESULTS AND DISCUSSIONS 


We compute various physical quantities using three 
small scale simulations, and compare them with the stan¬ 
dard results (STD). Disagreement with STD indicates 
the simulation strategy is incapable of capturing the cor¬ 
rect physics of infinite systems. 

(1) Average ion densities inside cavity 

We measure the average ion density in the simula¬ 
tion cavity as a function of the radial distance, shown in 
Fig. a If the multi-scale modeling is faithful, the result¬ 
ing average ion density must be independent of radius. 
In another word, non-uniformity of ion density indicates 
that the artifacts of cavity boundary have not been prop¬ 
erly cancelled. In Fig.^ we see that RP-I-GCMC simula¬ 
tion yields a much flatter density profile than RP-I-CMC. 
For the latter case, there is a systematic tendency that 
the ion density is higher near the center than near the 
boundary, with an overall variation of about one percent. 
This is caused by the artificial constraint of charge neu¬ 
trality inside the simulation cavity. 

Note that both RP-I-GCMC and RP-I-CMC show a 
thin boundary layer (with thickness comparable with ion 
diameter) where the ion densities vary rather abruptly. 
This is due to the short scale depletion effects of the hard 
wall boundary, and shows up actually in all simulations 
with hard wall boundary conditions. Even though cor¬ 
rection of these depletion effects are straightforward, it 
is more convenient just to exclude the boundary layer 
completely for data collection. 

This comparison clearly demonstrates that we should 


FIG. 5: Correlation potential acting on a positive ion, as a 
function of its position. The classical Debye-Hiickel theory 
works well in the low density regime. The substantial errors 
in RP-I-CMC arise to the artificial constraint of charge neu¬ 
trality inside the simulation cavity. The small difference be¬ 
tween DH and standard simulation results is caused by charge 
correlation effects. 


use grand canonical, instead of canonical, ensemble when 
implementing the multi-scale simulation strategy. 

(2) Correlation potential 

We measure the mean potential acting on an ion due 
to all other ions, which is usually called the correlation 
potential. Because we are in the low density regime, the 
Debye-Hiickel theory is applicable. Indeed, as shown in 
Fig. both RP-I-GCMC and STD agree well with pre¬ 
diction of DH. By contrast, RP-I-CMC shows substantial 
deviation from the other results. This is again due to 
artihcial constraint of charge neutrality inside the cavity. 

This correlation potential is not measurable in 
PBC-I-GCMC simulation, because the summation over 
potential due to all periodic images is divergent. 

(3) Screening charge density around a test ion 


We hx an ion inside the simulation cavity and measure 
the screening charge density around, as a function of the 
distance to the hxed ion. This charge density can be ob¬ 
tained using data of pair correlation functions inside the 
whole simulation cavity. The results have already been 
presented in Fig 1(a) Again, RP-I-GCMC agrees with 


STD up to high precision, whilst both PBC-I-GCMC and 
RP-I-CMC yield results that are substantially different. 

Theoretically, LPB predicts that the charge density has 
the form of screened Coulomb: 


Pg(r) =-eV^()>(r) = +qS{r). (55) 


Hence if we plot rpq{r) v.s. r in log scale, we should ob¬ 
tain a straight-line, with slope given by the inverse Debye 
length. As shown in Fig |l(a)[ this agrees remarkably well 
with the simulation results using RP-I-GCMC. 
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FIG. 6: The screening length as determined from large scale 
simulation and from multi-scale GCMC simulation. The solid 
curve is the prediction of PB theory. 


FIG. 7: A test ion is fixed at r inside cavity, and the total 
mean charge in the cavity is plotted as a function of r. Both 
STD and RP-I-GCMC agree well with LPB, Eq. (56l. 


We also plot the the analytic result of LPB with PBC, 
Eq. § in Fig |l(a)| as red diamonds (PBC-I-LPB). This 
agrees very well with Ewald method in the far field. Since 
Eq. ([^ differs from Eq. ( [S^ only by the periodic im¬ 
ages and their screening cloud, the latter can be identi¬ 
fied a main source of errors in simulations using periodic 
boundary conditions. 

The charge densities obtained using both RP-I-GCMC 
and STD can be used to determine the screening length as 
a function of ion density. The results are shown in Fig. 
Also shown there is the prediction of screening length by 
LPB, Eq. (36cI, which agrees rather well with simulation 
results. This shows that LPB is indeed applicable for the 
ion densities studied in this work. 


(4) Internal energy per ion 

We use different methods to compute the internal en¬ 
ergy of the system per ion. The results are shown in 
Fig |l(b)[ The total internal energy of the system can 
be easily calculated in terms of the correlation poten¬ 
tial acting on each ion: E = Again our 

multi-scale GCMC simulation yields the same results as 
the large scale simulation. Both multi-scale CMC and 
PBC-fGCMC yield incorrect results for the internal en¬ 
ergy. In the same figure, the solid line (DH) is the pre¬ 
diction of Debye-Hiickel theory in open space (the first 
term in RHS of Eq. which agrees with STD and 

RP-I-GCMC very well. The dashed line is the prediction 
of Debye-Hiickel theory with PBC as well as leading or¬ 
der correction of hardcore taken into account (Eq. (lOl). 
It appears that Eq. (101 takes into account most of the 
artifacts in PBC-fGCMC simulation if the density is suf¬ 
ficiently low. As the density increases, Eq. (101 deviates 
from PBC-fGCMC more significantly. 


(5) Total charge inside the simulation cavity 

As a stringent test of the faithfulness of our multi-scale 
method, we fix a positive ion at a distance r from the 
origin, and measure the total charge inside the cavity. 


The geometry is illustrated in Fig. The total mean 
charge inside the cavity is nonzero, because screening can 
not be perfect in a finite volume. In the framework of 
LPB, the total charge density is given by Eq. (551. Hence 
the total charge inside the cavity can be easily calculated: 


AQ(r) 


q- 


qK 

47r 


o-K|x-r| 


-d^x 


'\x\<R 


X - r 


qe -I-1) sinh(«:r) 

Kr 


(56) 


This result, shown as the solid solid curve in Fig. [V] 
agrees remarkably well with STD and RP-I-GCMC. 
Again, this indicates that 1) LPB works well in the den¬ 
sity regime under study; 2) our RP-I-GCMC simulation 
method faithfully captures the physics of an infinite elec¬ 
trolyte. By contrast, for RP-I-GCMC, this total charge is 
identically zero, due to the condition of charge neutrality. 


(6) Ion densities near a charged colloid 

Finally to test the applicability of our multi-scale sim¬ 
ulation strategy for inhomogeneous systems, we insert a 
uniformly charged colloid with charge Qq = 20e and ra¬ 
dius i?o = lOA at the center of the simulation cavity, and 
measure the local potential ^(r, q) acting on an ion fixed 
at r > i?o, due to both the charged colloid as well as all 
other ions. We set the dielectric constant of the colloid 
to be the same as that of the solvent, so that there is no 
image charge effect. Note that (p{r,q) also depends on 
the ion diameter d. We do not show this dependence in 
order to avoid cluttered notations. 

Note that the potential (p{r, q) is different from the av¬ 
erage potential (()(r) in the absence of the test ion. We 
shall call the former the conditional mean potential, and 
the latter the unconditional mean potential. Both quan¬ 
tities can be directly measured using simulation data. 
The difference between them arises because the test ion 
influences the distribution of other ions. The difference 
between (j){r, q) and (j){r) encodes essential information 
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FIG. 8: Potential and ion density profile near a charged spherical colloid, (a) <j>o{r), and (b) 0i(r) —(^i(oo). These potentials are 
defined in Eq. (57l. (c) Net charge density around the charged colloid. The decay length scales of (a) and (c) are approximately 
the same, and are about twice of that of (b), in agreement with theoretical predictions. 


about ionic correlations in this inhomogeneous system. 
We shall numerically compute q) and clarify its re¬ 
lation with (j){r). This demonstrates that our multi-scale 
reaction potential GCMC is capable of capturing the 
ionic correlations faithfully in an inhomogeneous system. 

The conditional mean potential 4>(r^ q) can be ex¬ 
panded in terms of the charge q: 

Hr,q) = Mr)+q(l>i{r)+ 0{q^). (57) 

It is important to realize that 4>o{r) is not the same as 
the unconditional mean potential 4>{r). Their difference 
is caused by the hardcore of the test ion. Nevertheless, in 
the far field regime where LPB is applicable, we can show 
that ipoir) and 4>{r) are proportional to each other, and 
both assume the form of screened Coulomb potential: 


Hence (j)o{r), (j)i{r) can be co mpu ted di rectly using sim¬ 
ulation data. In Fig. |8(a)| and |8(b)[ we plot repair) 
and r[(piir) — (()i(oo)] as a function of r in log scale. 
The far field asymptotics of both functions indeed be¬ 
have as screened Coulomb, with their characteristic de¬ 
cay lengths differ by a factor of two. This demonstrates 
that in our multi-scale CCMC simulation, the influences 
of the boundary of simulation cavity are properly can¬ 
celled, so that the simulation data can be used to analyze 
subtle ionic correlations. Also the fact that log (r(()i(r)) 
scales linearly with r in the far field indicates that the 
slow function x(r) is indeed changing very slowly. 


The potential of mean force (PMF) of a test ion (with 
positive charge q) can be obtained from the local poten¬ 
tial Eq. (57) via the method of Debye charging: 


p-Kr 

Mr) Air) -■ (58) 

r 

The potential (piir) is not completely known. In refer¬ 
ence [81] , it was analytically calculated near a uniformly 
charged fiat surface, by perturbing around the nonlinear 
PB theory. It was found that (piir) has the form 

_ 2 /^ 7 ^ 

(piir) = (piioo) + fir) x - -, (59) 

r 


where /(r) is a slowly varying function (comparing with 
exponential), and (^i(oo) is the reaction potential in the 
bulk. This result is expected to be valid in the present 
spherical geometry. 

Now using Eq. (57) we can express (pair), (piir) as lin¬ 
ear combinations of (pir, +q) and epir, —q): 


Mr) 

Mr) 


M, +q) + (pir, - 9 ) ) 

M, +q) - M, -q) 


(60) 

(61) 


U (r, q) 


(pin q) - 4’M,q) 


qMr) + Oie-^M 


dq + C/(r, 0). 


(62) 


The constant of integration U (r, 0) is the free energy cost 
of bringing a neutral particle from the bulk to r. In 
a symmetric electrolyte, U (r, 0) decays twice as fast as 
(pair), and therefore makes no contribution to the lead¬ 
ing order far field asymptotics of Uir,q). Now the aver¬ 
age density p+(r) of positive ions is related to the PME 
U (r, q) via the following well known Cibbs distribution. 
Hence the total charge density is 


pir) = qpa 


^-PU{r,q) _ ^PU{r-q) 


-2q^PMr)- (63) 


That is, the deviation of ion density from its bulk value 
also behaves as a screened Coulomb in the far field. In 
Eig. 8(c) we plot the corresponding simulation results and 
show that indeed p(r) obeys the law of screened Coulomb, 
with the same characteristic decay length 1 /k. 
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VI. CONCLUSION AND 
ACKNOWLEDGEMENT 

We have demonstrated the artifacts due to periodic 
boundary condition in the simulation of electrolytes. 
These artifacts are caused by the periodic image charges 
and the constraint of charge neutrality inside simulation 
domain, both unphysical from the perspective of mim¬ 
icking an infinite system. One can always avoid these ar¬ 
tifacts by simulating a much larger system and throwing 
away data from the region near the boundary. This leads 
to waste of computational resource, but is frequently 
done in simulations. Our multi-scale reaction potential 
GCMC provides a more efficient alternative. Combining 
reaction potential modeling and grand canonical Monte 
Carlo, this method cancel most of the electrostatic arti¬ 


facts introduced by the boundary of simulation data, and 
hence faithfully capture the physics of an infinite system 
even in a very small scale simulation (with linear size 
about three Debye lengths). If we combine this multi¬ 
scale CCMC method with appropriate fast algorithm for 
evaluation of electrostatic energy (such as the oct-tree al¬ 
gorithm), we arrive at a method for large scale MC sim¬ 
ulation of electrolytes, that will constitute a competitive 
alternative to the current methods based on molecular 
dynamics. This will be explored in a future publication. 
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